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We show how the band structure and beam dynamics of non-Hermitian PT- 
symmetric sinusoidal optical lattices can be approached from the point of view of 
the equivalent Hermitian problem, obtained by an analytic continuation in the trans- 
verse spatial variable x. In this latter problem the eigenvalue equation reduces to 
the Mathieu equation, whose eigenfunctions and properties have been well studied. 
That being the case, the beam propagation, which parallels the time-development of 
the wave-function in quantum mechanics, can be calculated using the equivalent of 
the method of stationary states. We also discuss a model potential that interpolates 
between a sinusoidal and periodic square well potential, showing that some of the 
striking properties of the sinusoidal potential, in particular birefringence, become 
much less prominent as one goes away from the sinusoidal case. 

PACS numbers: 42.25.Bs, 02.30.Gp, 11.30.Er, 42.82.Et 

I. INTRODUCTION 

The recent surge of interest in quantum Hamiltonians which are not Hermitian but which 
nonetheless possess a completely real energy spectrum, due to an unbroken PT symmetry, 
stems from the pioneering paper of Bender and Boettcher[l|], in which they showed, by 
numerical and asymptotic analysis, that the entire class of Hamiltonians 

H = p 2 - (ix) N (1) 

had that property for N > 2. Apart from the trivial case N = 2, the simplest example is 
for N = 3, where H = p 2 + ix 3 . 

Since that initial paper, there has been intensive investigation into the properties of 
Hamiltonians of this kind, whose progress can be followed in the reviews by Bender j2] 
and Mostafazadehj^]. We restrict ourselves here to those features that form the essential 
background to the present paper. 

For a viable framework of quantum mechanics one needs not only a real spectrum but also 
a probabilistic interpretation. In standard quantum mechanics that is provided by matrix el- 
ements of the type J ip*Ax, orthogonality of eigenfunctions J dx iplip2 = 0, and the probabil- 
ity density ip*ip. In PT-symmetric quantum mechanics it is found instead that orthogonality 
of eigenfunctions takes the nonlocal form J dx ipl(—x)ip2(x) = 0, that is, J dx (ipi)pTip2 = 0. 
In the context of quantum mechanics this is a problem because the metric involved in the 
corresponding normalization integral J dx (ip)pTip is not positive-definite, and we do not in 
the first instance have a proper probabilistic interpretation. However, it was subsequently 
found |I| that another metric could be constructed, with the help of a grading operator C, 
which preserved orthogonality and gave a positive normalization integral J dx (ip)cPTip- in 
contrast to standard quantum mechanics, this metric is not universal, but is dynamically 
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determined by the particular Hamiltonian in question. The calculation of this metric is 
usually extremely difficult, and in most cases can only be performed approximately, either 
through a set of algebraic relations or through the use of Moyal brackets 0]. 

A more general framework, of which CPT-symmetry is a special case, was developed by 
Mostafazadeh[3]. A Hamiltonian H is said to be quasi- Hermit ian if it can be related to a 
Hermitian Hamiltonian h by a similarity transformation: 

H = p~ l hp, (2) 

where p is a positive-definite Hermitian operator. From this we immediately see that 

H ] = r]Hr]-\ (3) 

where r\ = p 2 . The connection with the CPT formulation is that r] can be identified as 
when CP is written^ in the exponential form CP = e®. Accordingly p can be written as 



p = e 2^. 

The corresponding action of the similarity transformation on states is just 

\il>) = e&\<p), (4) 

where is a state of the non-Hermitian system governed by H and \<p) is the corresponding 
state in the Hermitian system governed by h. Then in the H system the positive-definite 
metric given by i] corresponds to the standard quantum-mechanical metric in the h system. 
Thus, [8|] 

<V>i|e-^2> = <Vi|e~* (e~*°e^) e -3G|^) = (5) 

A sur pris ing recent development has been the application of these ideas to classical 
optics [9|]-[l7j|- That such a transfer is possible is due to the fact that under certain ap- 
proximations the equation of propagation of electromagnetic waves reduces to the paraxial 
wave equation, which has the same form as the Schrodinger equation, but with different 
roles for the objects appearing there. The equation takes the form 

# = + rwU (6) 



dz \dx 

where now ip(x,z) represents the envelope function of the amplitude of the electric field, 
where z is a scaled propagation distance, and V(x) is the optical potential, proportional to 
the variation in the refractive index of the material through which the wave is passing. A 
complex V corresponds to a complex refractive index, whose imaginary part represents either 
loss or gain. In principle the loss and gain regions can be carefully configured so that V is 
PT symmetric, that is V*{x) = V(—x). There is also a non-linear version of this equation, 
arising from sufficiently intense beams, where there is an additional term proportional to 

Among many rece nt p apers we may mention linear (l3. and non-line ar|l7l| two-channel 



problems, and linear 111. 1 14. 1 151] and non- linear [l0| optical lattices. Ref. where, apart 
from an overall additive constant, the periodic optical potential was taken to be of the 
form V = |t4(cos2x + 2iVosin2a;), is of particular interest for the present paper, since it 
is a potential for which the equivalent Hermitian Hamiltonian can readily be constructed. 
Figures of the propagation profiles have been given, both below and above the threshold for 
PT-symmetry breaking at Vq = |, showing unusual features, such as non- reciprocity, power 
oscillations and bifurcation. In what follows we attempt to cast light on these phenomena 
from the point of view of the equivalent Hermitian system. 
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II. EQUIVALENT HERMITIAN HAMILTON! AN 



For the potential used in Ref. [11] . the analogue Schrodinger equation takes the form 

- ip" - - A(cos 2x + 2iV sin 2x)ip = -(3ip (7) 

for an eigenstate of H, with eigenvalue and z-dependence ip oc e~ ll3z . Below the threshold 
for PT-symmetry breaking, Vq < |, the real and imaginary parts of the potential can be 
combined into a cosine of complex argument, according to[l8[: 

cos 2x + 2iVq sin 2x = 7(1 -4V 2 )cos(2a;-^), 

where 9 = arctanh(2Vo). Thus, the non-Hermitian Hamiltonian 

H = p 2 — - A(cos 2x + 2iV sin 2x) (8) 

can be converted into the equivalent Hermitian Hamiltonian 

h = p 2 - - 4K 2 ) cos 2x (9) 

by the complex shift x — > x + \iO. This can be implemented by the similarity transformation 
of Eq. (J2]) , namely 

h = e -^ Q He^ Q (10) 

with Q = 0f> = —i9d/dx, which ensures that the spectra of the two Hamiltonians are 
identical. 

In the symmetry-broken case Vo > n> the corresponding identity is instead 
cos2x + 2iVosin2x = iy/{WQ - 1) sin(2x - <), 

where ( = arccoth(2Vo). However, in this case we have not gained a great deal from the 
similarity transformation, since the equivalent Hamiltonian h is itself non-Hermitian. 

Finally, at the critical value V = |, the equivalent potential vanishes altogether, so that 
the equivalent theory is simply a free theory, with spectrum (3 = —k 2 , as has been noted 
by Longhiflil]. among others. Part of this spectrum can be observed, in the reduced zone 
scheme in Fig. 1(b) of Ref. 0. The transformation in this case is a singular one, with 
6 — )• oo, so the methods used below can not be implemented for this limiting case. 



Band Structure 



In what follows we shall choose A = 4, the value taken in Ref. [Uj. Then the analogue 
Schrodinger equation for h for V < | is the Mathieu equation[l9|: 



ip" + (a — 2q cos 2x)ip = 0, 



(11) 



with q = — a/(1 — 4Vq) and a = —ft. In general terms the energy levels can be found by 
the Floquet method, whereby we take two independent solutions U\[x) and u 2 (x) with the 
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respective initial conditions ui(0) = 1, ^(0) = and w 2 (0) — 0, u' 2 (0) = 1 and integrate up 
to the Brillouin zone boundary at x = it to form the discriminant 

D(P)= 1 -(u 1 (n) + u' 2 (7r)). (12) 

If \D\ < 1, there exists a periodic Bloch-Floquet solution of the form 

¥k{x) = c k ui{x) + d k u 2 (x), (13) 

satisfying 

<p k (x + 7r) = e ifc > fc (x), (14) 

where k — (1/vr) arccos-D. This procedure gives fcasa function of (3, a relation that has to 
be inverted to give the band structure (3 = (3(k). In the standard notation for the Mathieu 
equation, k((3) is called the characteristic exponent. The values of /3 where k is an integer r, 
i.e at the Brillouin zone boundaries, are termed characteristic values, and are of two types, 
a r or b r , depending on whether the Bloch wave-function is even or odd. 

In fact in Mathematica these functions have been extended to k non-integral (the func- 
tions MathieuA and MathieuB), effectively mapping out the whole band structure (3{k) 
without the need to go through the Floquet procedure explicitly. Using this method we 
show in Fig. 1 the band structure in both the reduced and extended zone schemes for the 
potential of Eq. Q, or Eq. ©, for V = 0.45. 




FIG. 1: Band structure for Vq = 0.45 in the reduced and extended zone schemes. In the interests 
of clarity the gaps at \k\ = 2 have been slightly exaggerated. 

For Vo > 0.5 the analogue Schrodinger equation for h is 

ip" + (a + V(4V 2 - 1) sin 2x)^ = 0, (15) 

which, by a shift of x —> s — ir/2, again becomes the Mathieu equation, but with q pure 
imaginary. The functions a(k) and b(k) are still defined in Mathematica, apart from some 
minor glitches, and can again be used to map out the band structure. In this case, the 
characteristic feature, first observed in jlof, is that for real energies, which one is naturally 
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led to consider in solid-state physics, the bands no longer extend to the Brillouin-Zone 
boundary, but fold back on themselves. In the optical context, however, complex values of 
a are meaningful, corresponding simply to exponential growth or decay of the beam with z. 
The band structure for Vq = 0.7, derived using the same method, is shown in Fig. 2. 
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FIG. 2: Band structure (real part) for Vq = 0.7 in the reduced zone scheme. 



B. Bloch Wave- Functions 

In the present case the Floquet functions u\ and u 2 are precisely the even and odd Mathieu 
functions ce(a,q, x) and se(a, q,x) respectively, up to a normalization factor. For a given 
value of k, we know a and hence can determine the value of the ratio c k jd k in the equation 
( II 3p for the Bloch wave-function. Away from the Brillouin zone boundaries neither c k nor c4 
vanishes, but precisely at those boundaries one or other is zero, making the solution purely 
symmetric or antisymmetric. Thus, for example, at k — 1 the wave-function corresponding 
to the lowest band is symmetric, whereas the wave-function corresponding to the next band 
is antisymmetric. 

The Bloch wave-functions can be individually normalized in the extended zone scheme 
according to 

Vk {x)\ 2 dx=l (16) 

For k k! the orthogonality arises from the different periodicities of (fk{x) and ip^^x). If 
we use periodic boundary conditions in — Nit < x < Nn, so that k — » k r = r/N, 

[ <pt(x)<p ka (x) dx = e~^ A f ^Jp^ \ I ^ {x) Vks {x)dx, (17) 

where A = k r — k s = (r — s)/N. The second factor gives the orthogonality for r ^ s mod 2, 
while the remaining integral gives the orthogonality at the BZ boundaries. 
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III. METHOD OF STATIONARY STATES 

In quantum mechanics a standard method of implementing time development is the 
method of stationary states. That is, the initial wave-function ip(x,t = 0) is expanded as a 
superposition of orthonormalized energy eigenstates ipi(x): 

(p(x, t = 0) = Cj(pi(x), (18) 

i 

with 

and then 

<p(x,t) = J2 c Wi( x ) e ~ Eit - ( 20 ) 

In the optical problem exactly the same method can be applied, with z taking over the 
role of t, and the eigenstates being the Bloch wave-functions. We first apply this method 
to the Hermitian problem of Eq. fl9]) and then show how it can be adapted to give the 
^-development for Eq. (|7]). 



A. Propagation in Hermitian Case 

The initial envelope <p(x, z = 0) = g(x) is to be expanded in terms of the (p k r (x), according 



to 



g( x ) = ^2c r ip kr (x), 



(21) 



with the coefficients c r given by 



Nn 



Vl r (x)g{x)dx 



(22) 



-JVtt 



Using the translational property of the Bloch wave-functions we can reduce the integration 
range to the standard cell to tt: 



vll x ) G ( x ) dx 



where G{x) = J2q=-N e mqkr g{ x + ?tt), and then (p(x, z) is given by 



(23) 



E 



c r (p kr (x)e 



-ia(k r ) 



(24) 



For definiteness let us take q(x ) to be a broad Gaussian, g(x) = e ( x '/™) 2 ; with w = 6n 
, a function used in Refs. [ll|,[14j. In this case the absolute values of the coefficients are 
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FIG. 3: Absolute value of the coefficients c(k) for the Hermitian case. 

given in Fig. 3, which reveals that they fall off rapidly with \k\, and are essentially negligible 
for \k\ > 3. They are concentrated around the even integers, reflecting the slowly- varying 
nature of g. Note that by convention the wave-function is taken to be cek(x) at positive 
integers k and se^(x) at negative integers. 

When the sum of Eq. (T241 is performed, the development of the intensity, shown in Fig. 4, 
shows no surprises: the beam, representing a Gaussian wave-front at normal incidence, 
essentially propagates straight ahead, with a small amount of lateral spreading. 




FIG. 4: Intensity pattern in the Hermitian case. 
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B. Propagation in non-Hermitian Case 

The situation, however, is very different in the non-Hermitian case. We now need to fit 
the initial Gaussian g(x) to the transformed wave-functions ijjk(x) rather than to the (fk{x). 
That is, 

9( x ) = ^2d r ip kr (x). (25) 

r 

But from Eq. (J3J), which in terms of wave-functions reads ipk r ( x ) — ^Pk r { x — this can be 
recast as 

g(x+^i6) = ^2d r ip kr {x), (26) 

r 

and then ip(x, z) = <p(x — \i6, z). 

In fact, for the parameters taken, the coefficients d r differ very little from the c r . However, 
because we are plotting \ip(x, z) | 2 rather than \(p(x, z) | 2 , the intensity pattern is very different, 



showing the characteristic birefringence and power oscillations first noted in Ref . [11] . Figure 
5 is essentially identical to Fig. 2(a) in that paper. It turns out that the asymmetry is 
primarily due to the contribution of the Bloch functions with \k\ w2. 




FIG. 5: Intensity pattern in the non-Hermitian case. 



IV. DISCUSSION 

We have shown that it is a simple matter to derive the band structure of non-Hermitian 
sinusoidal potentials of the type occurring in Eq. (J7J) using the equivalent Hermitian Hamil- 
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tonian of Eq. and the corresponding Mathieu and related functions built in to Mathe- 
matica. The dynamics (z-development in optics) can also be implemented using the method 
of stationary states. Although in practical terms the z-development of the original equation 
is more efficiently found by direct numerical integration using the split operator method and 
fast Fourier transform, it is hoped that the stationary state method helps to elucidate the 
difference between the Hermitian and non-Hermitian situations. In particular it turned out 
for the parameters we used that the expansion coefficients did not differ significantly in the 
two cases: the difference between Figs. 4 and 5 was overwhelmingly due to the fact that 
in the Hermitian case one was plotting \<p(x, z)\ 2 , whereas in the non-Hermitian case the 
relevant function was instead the continuation \ip(x, z)\ 2 = \(p(x — \6, z)\ 2 . 

One may ask whether the birefringence shown in Fig. 5 is a general feature of PT- 
symmetric potentials, or whether there is something special about the sinusoidal potential. 
In particular, does the feature persist for a periodic square-well potential (which in practice 
would be much easier to construct)? The answer seems to be in the negative, and the 
transition from the sinusoidal case to the square-well case can be neatly studied by using 
Jacobi sn functions, whose m parameter allows one to interpolate between the two cases. 
That is, we replace 



V = ^A(cos 2x + 2iV sin 2x) 



(27) 



with 



W 



-A 



mi | K{m)^m \ + 2iVo sn ( — K{m),m 



(28) 



rescaling the argument so that the period remains tt. The imaginary part of the potential 
is shown in Fig. 6 for the two cases m = 0.8 and m = 0.999. 



Im(W) 



Im(W) 





FIG. 6: Imaginary part of the potential of Eq. (a): m = 0.8, (b): m = 0.999. 



For m = we reproduce Fig. 5, but by the time m reaches 0.8 the birefringence is much 
less prominent (see Fig. 7, upper panel), and for m = 0.999, at which point the potential 
is essentially a periodic square- well, it has more or less disappeared (Fig. 7, lower panel). 
Since in intermediate cases the wave-functions are not well-known functions, these figures 
were produced using the original method of direct integration of the differential equation for 
the z-development. 
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FIG. 7: Intensity pattern for W (Eq. ([25]) ). Upper panel: m = 0.8, lower panel m = 0.999. 

A few final remarks about the limiting case Vq — 1/2: for this particular value the 
eigenf unctions for the non-Hermitian V are again known functions, in fact modified Bessel 
functions Jfe(v2 e tx ), so that one might hope that the stationary state method could here 
be used directly for the V itself. However, as was pointed out by Longhi[14], this is not 
possible because of spectral singularities, or the non-completeness of the Bessel functions at 
the B-Z boundaries k = integer. This is another manifestation of the singular nature of the 
similarity transformation between the non-Hermitian and Hermitian problem in this case. 
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